PRO MAKE_BGPS_GALPOS_PLOT, SUBSET=subset
  
  IF n_elements(subset) EQ 0 THEN subset = ''
  
  ;;==========================================================================
  ;; Routine to make comparison plots for galactic positions derived
  ;;   from the NEAR/FAR determination of the code.  This analysis
  ;;   routine is a heavily-modified version of INVESTIGATE_GALPOS.pro
  
  ;; Define the COMMON BLOCK
  COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
     lpstr,rho_hi,rho_h2,rho_star
  
  ;; Since we are using routines buried in BGPS_FFORE, compile now to
  ;; they are available.
  FORWARD_FUNCTION BGPS_FFORE_MODEL, BGPS_FFORE_LOOPIE
  RESOLVE_ROUTINE, 'BGPS_FFORE', /COMPILE_FULL_FILE, /EITHER
  
  do_Tfit = 0b
  
  restore,'irdc_dist_model/bgps_rb3'+subset+'.sav',/ver
  s = read_bgps_csv(/ver)
  n = n_elements(rb3)
  
  
  ;; Define arrays for TEMP_DEPENDENT_VARIABLES
  tau    = dblarr(n)
  f_data = dblarr(n)
  sig_f  = dblarr(n)
  
  ;; Set up some of the variables needed throughout the program
  chi = '!9' + string("143B) + '!X'
  gr_rho = '!9' + string("162B) + '!X'
  gr_tau = '!9' + string("164B) + '!X'
  propto = '!9' + string("265B) + '!X'
  ;;"
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF NOT exists THEN galactic_params 
  
  ;; Using !MW parameters -- Goes into COMMON block!
  R0 = !MW.R0
  d = dindgen(!MW.NBINS*5.)*!MW.BINSIZE + !MW.BINSTART
  corr = R0 / 8.5d3
  
  ;; Set up the plotting environment   
  myps,'./irdc_dist_model/analysis_plots/bgps_galpos_plot'+subset+'.eps',$
       ct=0,xsize=5,ysize=7
  rotc=0b
  
  ;; Do re-usable calculations
  GENERATE_LPSTR, ROTC=rotc, SUFF=suff, RTIT=rtit
  
  
  
  ;; There are currently 4 cases of the model
  nmods     = 4
  dfn       = strarr(nmods)
  thisnf    = bytarr(n)
  fitdist   = fltarr(n)
  ffore_vec = fltarr(n)
  
  cnames = ['YGB7','Lime Green','TG3','RED7']
  colors = cgColor(cnames)
  textcolor = cgColor('RYB8')
  print,colors
  thick=[6,5,4,3]
  
  
  
  tstring = 'FIT'
  
  message,'Dust temp = '+tstring,/inf
  
  ytit='Galactocentric Distance [kpc]'
  ytf = "(I0)"
  
  plotsym,0,2.0                 ;,/fill
  plot,[-10.5,1.5],[-10,13.5],/nodata,xtit='Galactocentric Distance [kpc]',$
       ytit=ytit,ytickformat=ytf,/xst,/yst,$
       tit=rtit,$
       xtickformat="(I0)",charsize=0.75,/isotropic
  
  chmin = fltarr(nmods)
  ndof  = fltarr(nmods)
  tdp   = fltarr(nmods)
  
  FOR j=0, nmods-1 DO BEGIN
     
     message,'Case #'+string(j+1,format="(I0)"),/inf
     
     dfn[j] = './irdc_dist_model/data/case'+$
              string(j+1,format="(I0)")+$
              '_irdc_TdFIT'+subset+suff+'.sav'
     
     restore,dfn[j]
     chmin[j] = bestnorm
     ndof[j] = dof
     
     message,'Minimum Chisq = '+string(bestnorm,format="(F0.2)"),/inf
     
     do_Tfit = 1b
     
     residj = BGPS_FFORE_LOOPIE( res, NEARFAR=thisnf , FFORE_VEC=ffore_vec)
     
     ;; Sort out near vs. far vectors, and derive fitdist
     far = thisnf & near = ~ thisnf
     fitdist = near * lpstr.kdn + far * lpstr.kdf
     
     ;; Make x,y position for plotting
     x = -fitdist/1.d3 * sin(rb3.l*!dtor)
     y = -R0/1.d3 + fitdist/1.d3 * cos(rb3.l*!dtor)
     FOR kk = 0L, n-1 DO $
        plots,x[kk],y[kk],color=colors[j],psym=8,symsize=rb3[kk].c_meas
     
  ENDFOR
  
  ;; Solar Circle and Tangent Circles
  ells = findgen(1001)/1.e3*180.
  scl = findgen(1001)/1.e3*360.
  xells = -(R0/1.d3 * cos(ells*!dtor)) *  sin(ells * !dtor)
  yells = + (R0/1.d3 * cos(ells*!dtor)) * cos(ells * !dtor) - R0/1.d3
  xsc = -R0/1.d3 * cos(scl * !dtor)
  ysc = -R0/1.d3 * sin(scl * !dtor)
  oplot,xells,yells,linestyle=1,color=cgColor('Olive'),thick=3
  oplot,xsc,ysc,linestyle=1,color=cgColor('Dark Slate Blue'),thick=3
  
  ;; l=30 line (through end of the Bar)
  plots,0,-R0/1.d3
  plots,-R0/1.d3*cos(30*!dtor),R0/1.d3*sin(30*!dtor),/cont,linestyle=2,thick=3
  ;; al_legend,position=[-R0/1.d3*cos(30*!dtor),R0/1.d3*sin(30*!dtor)],$
  ;;           box=0,charsize=0.75,['l = 30']
  
  ;; Sun and Galactic Center
  plots,0,0,psym=7,color=cgColor('Orchid'),symsize=1,thick=2
  plots,0,-R0/1.d3,psym=7,color=cgColor('Sea Green'),symsize=1,thick=2
  
  plotsym,0,0.5,/fill
  al_legend,/top,/left,textcolor=textcolor,$
            ['Case '+string(indgen(nmods)+1,format="(I2)")+' ('+chi+$
             '!u2!dred!n = '+string(chmin/float(ndof),format="(F0.2)")+')'],$
            box=0,charsize=0.75,psym=8,colors=cnames
  
  
  myps,/done
  
END
